using DifferentialEquations,WGLMakie, JSServe Page(exportable=true, offline=true)
Model 1
Construct and solve and ODE problem to calculate the change in the position $X$ of a weight on a spring over time. I have rearranged the equation from the book to model the position $X$ rather than the force acting on it.
function model1(du,u,p,t) M, K = p X,V = u du[1] = V du[2] = -K*X/M end M,K = 1.0,0.5 u0 = [0.5,1.0] p = [M,K] # construct an ODE problem prob = ODEProblem(model1,u0,(1.0,20.0),p) # solve using Tsitouras 5/4 Runge-Kutta method sol = solve(prob, Tsit5(), saveat = 0.01);
Calculate the change in position and velocity for any given current position and velocity. Confusingly enough, $X$ is in the $y$ direction. As it's a spring, it feels intuitve to model its displacement vervically.
points = vec(Point.([-1:0.1:1]...,[-1:0.1:1]'...)) f((x,y)) = (-K*y/M,x) ./ (sum(abs.((-K*y/M,x))) * 20) c((x,y)) = sum(abs.((-K*y/M,x)))
c (generic function with 1 method)
Put them together in a plot
fig = Figure(resolution = (1600,800)) lines(fig[1,1],sol.t,hcat(sol.u...)[1,:], axis = ( xlabel = L"t", ylabel = L"X" ), ) arrows(fig[1,2],points,f.(points) , color = c.(points), axis = ( xlabel = L"\frac{dx^{2}}{dt^{2}}", ylabel = L"X" ), ) fig[0, :] = Label(fig, L"\frac{dx^{2}}{dt^{2}} = -\frac{KX}{M}", textsize = 24) fig
Model 2
Construct and solve and ODE problem to calculate the change in $X$ over time. This time with a damping parameter $K_{1}$.
function model2(du,u,p,t) M, K, K₁ = p X,V = u du[1] = V du[2] = (-K*X- K₁*V)/M end K₁ = 0.5 u0 = [0.5,1.0] p = [M,K,K₁] prob = ODEProblem(model2,u0,(1.0,20.0),p) sol = solve(prob, Tsit5(), saveat = 0.01);
Calculate the change in position and velocity for any given current position and velocity
points = vec(Point.([-1:0.1:1]...,[-1:0.1:1]'...)) k = 2 m = 1 f2((x,y)) = ((-K*y- K₁*x)/M,x) ./ (sum(abs.(((-K*y- K₁*x)/M,x))) * 20) c2((x,y)) = sum(abs.(((-K*y- K₁*x)/M,x)))
c2 (generic function with 1 method)
Put them together in a plot with a slider to control the $K_{1}$ value.
App() do session::Session K₁_slider = Slider(0.0:0.01:3.0) fig = Figure(resolution = (1600,800)) t = sol.t # We need to map a few differnt elements to the slider. # I'm sure there's a more succinct way of doing this. # But the following code gets the job done. # map diferential equation solution to slider u = map(K₁_slider) do val Nprob = remake(prob, p = [M,K,val]) sol = solve(Nprob, Tsit5(), saveat = 0.01) return hcat(sol.u...)[1,:] end # map arrow directions to slider dirs = map(K₁_slider) do val K₁ = val f2((x,y)) = ((-K*y- K₁*x)/M,x) ./ (sum(abs.(((-K*y- K₁*x)/M,x))) * 20) return f2.(points) end # map arrow colours to slider cols = map(K₁_slider) do val K₁ = val c2((x,y)) = sum(abs.(((-K*y- K₁*x)/M,x))) return c2.(points) end lines(fig[1,1],t,u, axis = ( xlabel = L"t", ylabel = L"X" ), ) arrows(fig[1,2],points,dirs , color = cols, axis = ( xlabel = L"\frac{dx^{2}}{dt^{2}}", ylabel = L"X" ), ) fig[0, :] = Label(fig, L"\frac{dx^{2}}{dt^{2}} = \frac{-KX - K_{1}V}{M}", textsize = 24) fig slider = DOM.div("K₁: ", K₁_slider, K₁_slider.value) return JSServe.record_states(session, DOM.div(slider, fig)) end
K₁:
0.0